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Abstract 


A  class  of  iterative  methods  for  solving  the  system  of  partial  differential  equations 
stemming  from  hierarchic  models  of  the  heat  conduction  problem  in  orthotropic 
laminated  plates  is  analyzed.  The  analysis  leads  to  sharp  explicit  expressions  for  the  rate 
of  the  convergence  depending  on  the  shape  or  director  functions  used  in  the  derivations 
of  the  hierarchic  models. 
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INTRODUCTION 


The  problem  of  partial  differential  equations  on  thin  domains  is  a  basic  problem 
in  engineering.  We  mention,  as  an  example,  the  typical  problem  of  the  homogeneous  or 
laminated  plates  and  shells.  The  usual  theoiy  is  either  based  on  physical  models  (see, 
e.g.,  [1],  [2],  [3],  [4]),  the  asymptotic  models  (see  e.g.,  [5],  [6])  or  the  hierarchic  models 
(see,  e.g.,  [7],  [8],  [9],  [10])  which  reduce  the  original  3  dimensional  problem  to  a  coupled 
system  of  partial  differential  equations  in  two  dimensions.  I'he  main  idea  behind  the 
hierarchic  modelling  is  to  create  a  sequence  of  two  dimensional  problems  the  solutions  of 
which  converge  to  the  solution  of  the  original  three  dimensional  problem.  The  problem 
of  a*posteriori  error  estimation  for  this  dimension  reduction  was  addressed  in  [9]. 

The  main  computational  problem  arises  from  the  obvious  fact  that  the  solution  of 
the  higher  order,  hierarchical  formulation  is  much  more  computationally  expensive  than 
the  lowest  one.  Hence,  a  major  question  arises,  how  to  solve  the  two  dimensional  elliptic 
system  by  an  iteration  procedure.  The  method  proposed  here  falls  in  the  general 
category  of  subspace  correction  methods. 

The  convergence  of  (parallel  and  successive)  subspace  correction  methods  under 
fairly  general  conditions  on  the  space  splittings  is  standard  by  now  [11].  The  most  widely 
used  splittings  are  either  related  to  domain  decompositions  or  to  h-multilevel  methods. 
Frequently,  however,  sharp  estimates  of  the  convergence  factors  cannot  be  obtained  for 
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situations  of  engineering  interest. 

In  the  present  paper,  the  function  space  used  to  derive  the  hierarchical  plate 
model  is  split  with  respect  to  the  (spectral)  order  of  the  transverse  shape  functions  used 
in  the  derivation  of  the  hierarchical  plate  model.  Further,  a  Fourier  analysis  of  the 
resulting  subspace  correction  method  is  performed  and  explicit  and  sharp  expressions  for 
the  convergence  factors  of  the  subspace  correction  method  in  terms  of  the  transverse 
shape  functions  "directors"  are  obtained.  This  allows  in  particular  to  (computationally) 
optimize  the  shape  functions  used  in  the  dimensional  reduction  with  respect  to  the 
convergence  rate  of  the  subspace  correction  method. 

The  outline  of  this  pf^r  is  as  follows:  in  Section  2  we  present  the  heat 
conduction  problem  and  the  hierarchic  models.  In  Section  3  we  introduce  the 
hierarchical  subspace  correction  methods.  Section  4  introduces  the  tool  of  Fourier 
analysis  and  contains  the  main  results.  Section  5  discusses  in  more  details  the  Fourier 
analysis  of  the  subspace  corrections  in  dependence  on  the  basis  used  in  the  modelling. 
Section  6  finally  presents  numerical  experiments. 

2.  THE  BOUNDARY  VALUE  PROBLEM  AND  THE  HIERARCHY  OF  MODELS 
For  Oj  G  R“,  we  define  the  (hyper)  rectangle 

(2-1)  :«  (x  €  R"|  |x.|  <  a,/2,  1  £  i  ^  n}. 

We  will  be  mainly  interested  in  n  =  1  or  2,  but  our  analysis  is  for  convenience  for  any 
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n  €  N.  With  and  0  <  d  <  1  we  associate  the  domain 
(2  2)  0  fc,  X  (-d,d) 

of  thickness  t  =  2d  and  the  "edge" 

(2.3)  r  »»  7  X  (-d,d),  7  « 

Points  in  Q  shall  be  denoted  by  (x,y)  where  x  €  «^,  |y|  <  d.  We  also  define  the 
faces 

(2.4)  .  {(x,y)|x  €  y  *  ±d}. 

In  Q  we  consider  the  boundary  value  problem 

Lu  «  0  in  Q 

(2.5)  u  -  0  on  r, 

Du  «  f*  on  R^, 

where  D  denotes  the  conormal  derivative  and  L  is  defined  (in  the  sense  of 


distributions)  by 


L«  .  -ifafil  «!ll  .bU 

dy  [  [dj  dy j  [d 


where  V,  denotes  the  gradient  with  respect  to  x  €  R*  and  a,b  €  L*(-l,l)  satisfy 


(2.6)  0  <  a  ^  a(2),  0  <  b  s  b(z) 

and  C  =  CT  isan  nxn  matrix  for  which 

(2.7)  0  <  c  1(1*  S  rC(,  V  0  €  R-. 


For  the  weak  formulation  of  (2J)  define 
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H‘(0,r)  «  H‘(0)n{ulu-0  on  T). 
Then  (2 J)  becomes:  Find  u  e  H‘(Q,r)  such  that 


(2.8) 

B(u,v) 

-  F(v),  w  €  H‘(n,r), 

where 

(2.9) 

B(u,v)  -  1^1  a  1 

(2.10) 

F(v)  -  J 

{r(x)v(x,d)  +r(x)v(z,-d)}dx. 

Under  the  assumptions  (2.6),  (2.7)  the  problem  (2.8)  has  a  unique  weak  solution 
u  €  H’(Q,r)  provided  that,  for  exaiT'ple  f*,f*  €  V(taJ  (this  could  be  weakened). 
For  simplicity  we  assume  also  that 

a(z)  ■  a(-z),  b(2)  »  b(-z),  a.e.  z  E  (-1,1) 

and  that 

(2.11)  r  -  f-  -  f, 

which  ensures  u(x,y)  =  u(x,-y)  a.e.  in  Q. 

Hierarchical  models  of  (2.8)  are  approxiraations  of  (2.8)  by  elliptic  boundary  value 
problems  on  which  are  obtained  from  (2.8)  by  a  dimensional  reduction  method 
which  we  now  describe.  Let 
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u  .  £  X|(xV|(^),  X|(x)  e 

i-0  O 

with  linearly  independent  director  functions  ^  and  such  that  span 

*  span{\^j}jl4  where  ^j(z)  >  \^j(-z)  and  the  functions  are  linearly 
independent  and  recursively  defined  by 
(2-13)  j'^a(z)^Jv^d2  «  0, 

(2.14)  j  *^a(z)^(v'dz  +  |*^b(z)^jvdz  ■  v(l)  ♦  v(-l), 

(2.15)  J' a(z)^(v'dz  |U(z)^j.,vdz  -  0,  j  -  2^,...  vv6H'(-l.l). 

This  selection  of  the  functions  ensures  certain  optimality  properties  of  the 

hierarchical  models  [7].  Moreover,  it  was  shown  in  [7]  that 

(216)  is  dense  in  H‘(-l.l)  fl  -  ^-*)> 

Remark  2.1.  From  (2.13)  it  follows  that  is  constant. 

The  sequence  of  Neumann  problems  (2.14),  (2.15)  defines  the  functions 
uniquely  since  the  constant  for  is  uniquely  determined  by  the  compatibility  condition 
ensuring  the  existence  of  □ 

Obviously,  the  space  S(q)  in  (2.12)  is  a  closed,  linear  subspace  of  H’(Q,r)  and 
hence  the  dimensionally  reduced  problem  is  :  find  u^  €  S(q)  such  that 


(2.12) 


S(q) 
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(2.17) 


B(u„v)  =  F(v),  vv  €  S(q). 


This  problem  admits  a  unique  solution  u^.  The  problem  (2.17)  results  in  the  following 
boundary  value  problem  for  X  :=  (Xj,X,,...,X )’’’  6  (Hi(wJ]’*‘: 


(2.18) 


L(D^)X  :=  -d*AP(D,)X  BX  =  2df(x)^(l)  in  «, 


X  =  0  on  du_ 


where  ^z)  :=  {<9^,^,,...,^ and 


(2.19) 


A  =  J ’^b(z)^^‘*'dz,  B  =  j’^a(z)v>V' dz 


and  the  differential  operator  P(DJ  is  given  by 


(2.20) 


P(DJ  =  V,.CV,. 


Remark  2.2  More  generally  than  (2.12),  S(q)  may  be  defined  to  admit  a  model  order 
q  which  varies  throughout  This  will  not  be  considered  here  and  we  refer  to  [13], 
[14]  for  details,  □ 


Remark  2.3.  Other  boundary  conditions  on  F  could  of  course  also  be  considered.  □ 


Remark  2.4.  In  general,  the  boundary  value  problems  (2.18)  must  be  solved 
approximately,  e.g.,  by  the  finite  element  method,  but  we  will  not  address  this 
approximation  explicitly  and  assume  that  all  boundary  value  problems  on  u,  are  solved 
exactly.  □ 
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Remark  2.5.  Define  V(q)  :=  span{^j}j|^  C  Then  obviously 

S(q)  «  hJ  (wJ  0  and  the  dimensionally  reduced  solution  u,  in  (2.17)  is 

independent  of  the  particular  basis  of  used  in  the  definition  (2.12)  of  S(q). 

□ 

3.  HIERARCHICAL  SUBSPACE  CORRECTION  METHOD  (HCM) 

Solving  the  problem  (2.18)  by  the  finite  element  method  (cf.  Remark  2.4),  a  large 
system  of  linear  equations  has  to  be  solved.  We  formulate  therefore  block  iterative 
techniques  for  the  solution  cf  these  systems,  each  block  corresponding  to  one  equation  in 
(2.18),  with  a  size  comparable  to  that  of  the  simplest  model  in  the  hierarchy  for  q  =  0. 

We  assume  here  that  the  equations  (2.18)  of  the  hierarchical  model  have  been 
discretized  with  high  accuracy.  This  allows  to  formulate  and  analyze  the  iterative  scheme 
in  the  semidiscrete  setting,  i.e.  under  the  assumption  that  (2.18)  is  solved  exactly. 

All  methods  considered  can  also  be  interpreted  as  subspace  correction  methods 
(see,  for  example,  [11])  based  on  a  suitable  splitting  of  S(q)  in  (2.12)  which  we  now 
introduce:  we  write 

(3.1)  S(q)  -  S„  S,  +  • . .  ^  S, 

where 

Sr'HJw® 

and 
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with  as  in  (2.13)-(2.1S).  Since  the  are  by  assumption  linearly  independent,  we 
have 

(3-2)  n  S„  =  {0},  m  ;«5  n,  0  s  m,  n  s  q. 

We  can  therefore  write 


and  (2.15)  takes  the  form 

(3-3)  £  B(Uf.vJ  =  F(vJ,  vv,  €  S„  k  =  0 . q, 

i-o 

s 

where  v  «  X)  €  S(q),  €  S^.  Evidently,  (3.3)  and  (2.17)  are  equivalent  weak 

klo 

formulations  of  the  reduced  problem  (2.18). 

Our  subspace  correction  algorithms  will  be  based  on  the  decomposition  (3.1).  We 
begin  by  formulating  a  successive  subspace  correction  algorithm  with  relaxation  which 
becomes,  upon  FE  discretization  of  (2.18)  in  the  x  variable,  a  block  SOR-method. 
Algorithm  3.1.  1.  Given  a  relaxation  parameter  0  <  0  <  2,  and  an  initial  approxima¬ 
tion  u®.  Set  n  =  0. 

2.  Repeat  until  convergence  based  on  a  given  tolerance: 

2.1.  u“  *  ^  U|"\  U;">  G  Sj 

j*o 

2.2.  For  k  =  0,...,q  solve  the  problem:  Find  €  S^  such  that 
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(3.4) 


B(ur‘\vJ  »  (1-0)b(u<”>,vJ 

-  B(urv.)  *  i;..  b(u»’4 

+  tfF(v^),  Vv^  G 


2.3.  Set 

u(-»  = 

i*o 

n  :=  n  +  1. 

2.4.  end. 

We  already  observed  in  Remark  2.5  that  there  is  considerable  freedom  in  the 
selection  of  the  basis  functions  which  does  not  affect  u^.  It  is  apparent,  however, 
that  the  space  decomposition  (3.1)  strongly  depends  on  the  basis  chosen  in  (2.12). 
This,  in  turn,  shows  that  the  convergence  properties  of  Algorithm  3.1  are  governed  1^ 
the  basis  Our  purpose  in  the  following  sections  is  a  quantitative  analysis  of  the 

convergence  rate  k  of  the  Algorithm  3.1  depeniuig  on  the  basis  and  the  data 
f(x)  of  the  problem.  By  convergence  rate  we  mean  the  smallest  number  i:  G  [0,1]  such 
that  for  any 

where  ||v||g^o,  -  (B(v,v))''^. 

The  convergence  rate  k  in  (3.5)  guarantees  a  reduction  o';  the  error  by  a  factor 
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K  per  step.  Frequently,  however,  the  asymptotic  convergence  rate 


1/k 

(3*6)  f,  *  sup 

Hm 

”  "  E(fl) 

fSLV.) 

k-*oo 

is  a  better  measure  for  the  actual  performance  of  the  HCM. 

Remark  3.1.  As  we  will  see  in  the  following  section,  we  always  have  p  s  *.  In  our 
experience  the  observed  rate  of  convergence  is  much  closer  to  p  than  k  except  for 
possibly  the  very  few  first  steps.  □ 

Remark  3.2.  Algorithm  3.1  could  be  interpreted  in  the  usual  fashion  as  a  multiplicative 
subspace  correction.  Since  it  is  based  on  the  hierarchy  of  models,  we  shall  refer  to  it  as 
Hierarchic  Subspace  Correction  Method  (HCM  for  short). 

Along  the  same  lines,  of  course,  also  additive  subspace  corrections  based  on  the 
splitting  (3.1)  could  be  considered.  □ 


4.  FOURIER  CONVERGENCE  ANALYSIS  OF  THE  HCM  -  I 

In  this  section  we  give  a  convergence  proof  for  the  HCM  using  Fourier  analysis  in 
the  case  of  problem  (2.5).  To  this  end  we  define  for  o  e  R"  the  (hyper)  rectangle 
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(4.1) 


w_  := 


a.  Sa- 

X  €  R"  I  -  <  X.  <  _' 

'  2  ‘  2 


A  function  w  on  R°  is  called  a  -periodic,  if 

w(x+o,e.)  =  w(x),  a.e.  x  G  R",  vi. 
The  vector  ^  corresponding  to  a  is  defined  by 


T  X  X 


To  present  the  weak  form  of  (2.5)  in  the  periodic  setting,  we  collect  first  a  few 
notions  on  periodic  function  spaces  (for  more  on  the  asymptotic  analysis  of  hierarchical 
plate  models  in  the  periodic  setting,  we  refer  to  [15]).  Throughout  the  present  section 
we  will  understand  all  Sobolev  spaces  as  spaces  of  complex  valued  functions  and  indicate 
by  a  bar  the  conjugate  complex  quantity.  By  H^(uJ  C  H'(a)J  we  denote  the  space 
of  all  2cr-periodic  functions  w  G  ^(ojj. 

If  w  G  then  let 

(4.2)  w(x)  =  (2x)-”  Y: 

■  6Z' 

where  the  Fourier  coefficients  w,  arc  given  by 

W.  -  (iSA-lSJ  f  w(x)ei«IP«>d;. 

and  [inj5]  denotes  (m,/S,,mj/3j,...,mj3^).  For  v(x)  =  (2x)"“  Y  we  also 

■€Z* 

have  Parseval’s  formula 
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(43) 


[  wvdx  =  (2t)-"(/3,/Sj.../3J-‘  £  w.v,. 
J  -ez* 


Fourier  transforms  of  vector  functions  are  taken  componentwise.  Analogously,  by 

we  denote  the  set  of  all  w(x,y)  G  which  are  a  -periodic  with  respect 

to  the  first  variable  in 


(x,y)  €  n,  :=  X  (-d,d). 
From  Parseval’s  equation  we  obtain  immediately 
Proposition  4.1.  For  all  u  G  per‘(ftj 


(4.4)  "LI 


dxdy 


.  (2T)-'(ff,...S.)-'  p  |'  j<|-'a(z)|«'.(2)p.db<2)[ii/jrqiiiffllu.(z)|>} 


The  expression  (4.4)  in  Proposition  4.1  is  a  norm  on  the  subspace 

>  H'CfljnlulLu^dy-o}. 

Further  define 

L=(ft.)  :•  n  {  “  1 1-  «<h  -  o}. 

Defining  Sp„(q)  ;=  ®  V(q)  and  Sp,,(q)  :=  H^(«.)  fl  we  can  write 

u  (x,y)  =  X(x)T 


^  ^  j  G  Sp^,(q)  in  the  form 


(4.5) 


u,  =  (2t)-"  x; 

■ez*  0 


where  Xm  ~  arc  (the  vectors  of)  Fourier  coefficients  of  order  m  of  the 
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(vector)  function  X(x).  Moreover 


(4.6) 


■ez* 


■67 


(AT)  (|niprqmPl)“ 

and  where  the  matrices  A  and  B  are  given  in  (2.19).  Then  the  symbol  of  the  differential 

operator  L(DJ  in  (2.18)  is 

(4.8)  q,J  ,iA  .  B 

Due  to  Proposition  4.1.  the  norm  defined  by  (4.6)  is  a  norm  on  the  subspace 

fl«r(<))- 

Remark 4.1.  If  <pg  »  const  and  =  0  for  j  S:  1,  we  find  that 

^  ^  ^per(^)  »  Xoc  '  0  in  (4..5). 

More  generally,  for  let  1  =  eV  for  some  e  €  S’*’.  Then 

u,  €  S'pgj.(q)  «  e’^x,  »  0  in  (4.5).  □ 

By  B^(»,*)  and  F^(»)  we  denote  the  forms  in  (2.9),  (2.10)  with  integrations 
over  0^  and  respectively  (and  using  v  instead  of  v). 

Proposition  4.2.  Let  f  €  L*(«J,  f(x)  =  (2t)‘”  f^e''**’^***^ .  Then 

■6Z‘ 

I  f(x)dx  =  0  ♦♦  f j  *  0  and  if 
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(4.9)  fo  =  0 

then  the  boundary  value  problem  :  Find  u  E  0 

(«»)  B.(ii.v)  .  F.(v),  vv  6  fl’pe/fij 

and  the  reduced  prtMem:  Find  E  S’pgj.Cq)  such  that 


(4.11)  B>,,v)  =  F^(v),  Vv  E  ?pej(q) 

both  admit  solutions  which  are  unique  up  to  constants.  □ 

We  shall  now  link  the  problem  (2^)  on  0  to  an  equivalent  one  on  6^.  This 
will  be  done  by  an  extension  operator  ^  defined  as  follows.  We  subdivide  into  2° 
hypercubes  and  identify  one  of  them  with  in  (2.1).  We  extend 
u  E  H‘(0,r)  to  fi  E  as  follows:  in  2"*’  octants  we  select  il  =  -u,  with 

the  argument  of  u  properly  shifted,  and  in  the  remaining  2°'*  -1  octants  we  select 
fi  -  u.  (We  note  that  u  =  0  on  dw,  X  (-d,d)  and  that  this  guarantees 
Q  E  H^(0j.)  In  the  same  way  we  extend  the  function  f.  Here,  of  course,  only 
f  E  L*(<J>J. 

Proposition  A3.  The  extension  :H'(0,r)  -►  H^(6j  is  continuous  and  isomemc,  in  the 
sense  that 


I"  HU  -  liri;, ,  - 

Proof:  Due  to  the  homogeneous  essential  boundaiy  conditions  satisfied  by  u  on  P  we 


observe  that  u  E  H‘(A^)  and  that  u  =  0  on  bCa^  x  (-d,d).  Hence  it  admits  an 
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a-periodic  extension  to  R"  x  {-d,d)  which  is  locally  in  H*.  This  implies  that 
u  E  The  properties  [  u  dxdy  »  0  and  (4.12)  are  obvious  from  the 

definition  of  ^  □ 

Remark  4.2.  We  have  actually  shown  the  extension  ^  to  be  continuous  fi’om 
onto  the  subspace  of  of  functions  the  traces  of  which  vanish  on  the 

hyperplanes  x.  =  a^/2,  i  =  1,2,. ..,n.  □ 

Remark  4.3.  The  space  ®  closed,  linear  subspace  of  Hje,(6,)  and  we 

define 

Vo)  >  5p<,,(q)  n  VV 

Then,  if  we  set  for  given  f  G  L*((i)  J 

F.(v)  ==  j.  ^ (x)(v(x,d)  +  v(x,-d))dx, 

(4.10)  and  (4.11)  admit  unique  solutions  0  G  A^(fij,  0^  G  Sp„(q).  Their 
restrictions  to  0  coincide  with  the  solutions  of  (2.8)  and  (2.17),  respectively.  □ 

Remark  4.2  is  the  basis  for  an  equivalent  formulation  of  the  HCM  Algorithm  3.1 
in  the  periodic  setting.  To  this  end  we  introduce,  analogous  to  (3.1),  the  subspace 
decomposition 

(4.13)  S^(q)  .  S,  .  S,  *  . . .  .  S, 

and  we  observe  that  ^:Sj  —  Sj  continuously.  Then  the  periodic  version  of  Algorithm 
3.1  is  defined  exactly  as  in  (3.4)  and,  due  to  Remark  4.3,  the  restrictions  of  the  resulting 
iterates  €  Sp^XQ)  to  are  those  produced  by  Algorithm  3.1,  i.e. 
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Vn. 


(4.14)  .  uf‘ 

From  this  equivalent  periodic  formulation  of  Algorithm  3.1,  however,  a  precise  estimate 

of  the  convergence  factors  p  and  k  in  (3J)  and  (3.6)  can  be  derived. 

Theorem  4.1.  Let  u^\x,y)  €  ^^(q)  denote  the  k-th  iterate  of  the  successive  subspace 
correction  edgoriihm  3,1.  Let  further  for  a  relaxation  parameter  0  <  $  <  2 

(415)  : -  (D  -6D  - «E«) 

be  the  SOR  iteration  matrix  corresponding  to  L(0  defined  in  (4.8)  (as  it  is  usual, 

D(i’)  and  E(f)  denote  the  diagonal  and  strictly  lower  triangular  parts,  respectively,  of 

UD).  Then 

i)  the  contraction  constant  k  in  (3.5)  admits  the  estimate 

(4.16)  s  sup(X(«,i?,d)) 

•6z; 

where  X(0,ii)  is  the  largest  eigenvalue  of 

M(tf,ij)‘*L(ij)M(fl,i7)x  -  X(0,t;)L(i})x, 
and 

z;  -  z-  n  e  a;„(«.)} 

ii)  the  asymptotic  convergence  rate  p  in  (3.6)  admits  the  estimate 

(4.17)  P  '  sup{r(M(^,ij,d))} 

■€z; 

where  r( . )  denotes  the  spectral  radius. 

Proof:  We  proceed  in  several  steps. 
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i)  We  write  \{x,y)  «  X(x)V|^ 
and  have  the  Fourier  expansionr 


,  -  X«(x)V 


X(x)  .  (2t)-J  E  X. 

-€z; 

X«(x)  -  (2t)<-">  E 


where 


i€z; 


Xm  *  (x2>.xi\.-.x£,)  . 


etc. 


Let  further 


r(x)  -  ft  -  (2t)-“  J; 

•€2; 


Then  the  Fourier  coefficients  x.  of  the  vector  function  x(x)  are  determined  from  the 
linear  system 

(4.18)  Mi|.d)x.  «  2df.^l),  m  €  ZS 

where 

V,  >  (lm^fqin«)>«. 

The  crucial  observation  is  now  that  the  successive  subspace  correction  algorithm 
3.1  becomes,  in  our  periodic  setting,  a  (frequent  dependent)  SOR-method  for  the 
solution  of  the  hermitean  linear  systems  (4.18),  i.e. 


M((?.i7.d)x.(k)+  2df.lCl(tf,,.d)^(l), 


where 


(4.19) 


lVl(fl,)?.d)  -  tf(D(,.d)  ..  eE(i,.d))-‘. 


In  particular,  no  Fourier  modes  that  are  absent  in  the  right  hand  side  f(x)  and  the 
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initial  approximation  X<*’(x)  will  be  excited  in  the  course  of  the  iteration. 

ii)  From  (4.10)  we  obtain  that 

-  E  (x.-xifMi.'fllit.-xJ’) 

■62; 

where  C,  -  d-‘(2T)'^/3,...^,)-‘. 

With  the  Fourier  coefiScients  of  the  error  X  - 

•“  Xa  ~  x2^  ®  ^ 

we  find 

|4,  -  Clirtj  •  c,  E  ivS'U 

where 

iY«n  .  vj'»4,,d)v? 

This  proves  (4.16), 

iii)  To  prove  (4.17),  we  first  establish  a  perturbation  result.  Recall  that 
L  =  ij*A  +  B,  A  and  B  as  in  (2.19)  and  (4.8).  Then 

D  -  +  D,,  E  «  7/%  ♦  E,, 

and  we  have  for  sufficiently  large  i) 
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M  -  ((d’D.+D,)  ♦  «(il’E,*E,))  '  {(!-»)(, ‘D.^D,)  -  «(VE.*E,)"} 

-  ((D,*»E.)  +  ,-=p,t»E,))"  •  {((1-9)D.-»E?)  .  ,-’((1-*)D,-«E!')} 

where  the  remainder  R  is  given  explicidy  by 


•  £  (  .,-‘p,.»E,)-'p,.«E,))'  •p.*«E.)-'|((l-«)D.-«Er)  *  ,--((l-9)D,-9E,")} 

j'* 

*  ,-‘p,*9E.)-'((l-«P,-9E,“) 

provided  that  v  >  Vo>  is  selected  sufficiently  large  to  ensure  the  convergence  of 
the  Neumann  series.  Hence  we  have  shown  that,  as  ->  «», 


(4.20) 


M(9,v)  -  M^(0)  *  R(e,i,) 


where  M^(0)  is  the  (ii  -  independent)  SOR  matrix  corresponding  to  A  and 
II  RH  *  0(i;'^)  as  ri ’■*  oa  in  any  fixed,  t;- independent  matrix  norm, 
iv)  We  prove  (4.17).  As  before,  we  find  for  k  €  N, 


u  -uji*  .  2-C„  J:  (v<«)  (MrU^.cl)M'T^r 


2-C,  S  ||Mnf<'>|*. 


where  ||Y||l«  :«  ||Y«LY|j^ 
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As  it  is  well  known,  for  eveiy  matrix  M  and  every  c  >  0  there  exists  a  vector 
norm  |||*|||m  on  O'*'  such  that  for  the  induced  l.u.b  matrix  norm  (which  we  again 
denote  by  |||  •  |||,^)  there  holds 

(^•21)  III  Mill  S  r(M)  -.-5/2,  6  >  0 

(the  norm  |||*|||m  depends  on  e). 

In  particular,  we  observe  that  there  exists  a  constant  c  ^  1 
independent  of  ij)  such  that 
(4.22a)  |||R(«,a)|||  s  c/,>,  , 

(4.22b)  r(MJ  -  c/if*  S  r(M)  S  r(MJ+  c/i»*. 

Now  we  fix  e  >  0  and  introduce  the  index  sets 

C.  =  {m  €  Z;|c/(i,id*)  <  €/4},  C,  -  2;\(J. 

with  the  constant  c  as  in  (4.22)  and  tiji  >  ri^.  Then  we  split 

.  2-c^  j  r  *  i:  j 

and  discuss  both  sums  separately.  For  m  €  Sf,  we  have 

-  (v«)''(M‘)"Ma.<i)Nm(r 

-  (CflNPHaLd’A  .  B)(M^vr 

S  aid’||A‘”M*V*||l  •  (**’)"(M‘)"BIiOTC. 

We  estimate  the  first  term  by 
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s  ,Ld>iA|AM)lll(M.-R)‘iiiL. lv.«n 

s  >id>A;M)a)nd,(A)(|||M,|||„_  »  II|R||Im.)“K“1I- 

wfaere  A,(6,d)  is  the  equivalence  constant  between  |  •  I,  and  III*  III  for 
sufficiently  large  k.  For  sufficiently  large  k,  the  second  term  can  be  estimated  by 

^  r(B)l|M''YJ*>||^Aj(c,d) 

S  nidh/4 .  r(B)  I  A‘«M^Y.<»>  ||lAi(A)A?(e,d) 

-  A,(€,d)i,id*||A‘«M^J®>|^. 

Then  we  continue  as  in  the  estimate  of  the  first  term.  By  (4.21),  this  yields  with  a  new 
constant  A(6,d)  the  estimate 

IhfvS’lL.  s  A(«,d)eond,(A)(t(M.).e/4.e/4)”lYriL„ 
and  hence  we  obtain  for  every  c>0,  0<d<l  and  k  € 

(4J3)  E  S  cond,(A)A(e.dXr(MJ  .  s/2)‘  E  K'||[, 

■ec.  Mec. 

We  consider  the  case  m  E  Since  is  a  bounded  set  for  fixed  e  and  d, 
we  estimate 
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E  IwvS’lL-  s  5;  |M‘||L4vri,^ 

met.  ■€«. 

(4^)  S  £  A.M)ll|M‘llii||'d‘t- 

'  ■'  ■€<[. 

£  (  sup  A.(4.d))5;  (i(M(«,,.d))  »  s)»|Y||i.  . 
mGC^  -es. 

From  (4.22b),  we  see  that 

r(M^(fl))  S  r(M(d,i,^d))  +6/2,  Vm  6  C„ 
so  we  can  combine  (4.23)  and  (4.24)  into 

|u,->>f’rE  -  E  (r(M(*,,.d))  *  t)» 

■€z; 

s  2'(C;,C.(.,d)(r-<)»  £  iyTIU 

■ez; 

-  C,M)(f  ».)»||u, -<>111, 

where 

(4.25)  r(«)  •  swp  r(M(e,i?,d)) 

mez; 

and  Q  depends  on  e  and  d,  but  is  independent  of  k.  Hence 

l/k 

lim  ^  (r+6)nm  (C(6,d))'* 

k-00  ||u,-u;'»||,J  k-eo' 

=  r  +  6. 

Since  e  >  0  and  0  <  d  <  1  were  aibitrary,  the  assertion  (4.17)  follows.  □ 

Remark  4.3.  In  (4.17)  we  have  taken  the  sup  over  all  ij„d  G  Z?.  Obviously  we  get 
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w 
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an  upper  bound  when  we  oefine 

5(^)  «  sup  {r(M(«,»f))}. 

D>0 

Theorem  4.1  has  several  important  consequences  which  we  discuss  next.  The 
crucial  observation  was  (4.16)  and  (4,17)  which  states  that  in  the  periodic  setting  the 
HCM  becomes  simply  a  (parameter  dependent)  SOR  method  for  the  solution  of  the 
linear  system  (4.10).  Since  for  every  m  f  0  the  matrix  L(i;.d)  is  symmetric  and 
positive  definite,  we  have  from  Kahan’s  theorem  (e.g.,  12,  Thm  83.5)  the 
Corollary  4.1.  There  holds,  for  every  d  €  (0,1]  and  m  € 

(^•26)  |a,-l|  ^  p  s  Jt  . 

Proof;  We  have  Kahan’s  lower  bound  |  oi  - 1 1  ^  p(M(0,i;))  for  all  it  ^  0  (see,  e.g. 

[12,  Chap.  8|)  and  hence  (4.26)  follows  from  «  ^  p(M(tf,ij)).  □ 

Corollary  4.2.  The  rate  k  of  convergence  of  the  HCM  is  bounded  from  above 

independentfy  of  d. 

Proof:  From  Theorem  4.1  we  have 

(«)*  S  sup  {X(®,^,d)} 

meZo" 

S  sup  {X(tf,i7)) 

where  X  is  as  in  (4.16).  *nCR 

The  result  follows  since  the  matrices  M(6,i})  and  L(i;)  are  independent  of  d.  □ 


The  SOR-iteration  applied  to  a  diagonal  matrix  L  converges  in  one  step. 
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Consequently,  if  we  can  give  a  basis  <p  that  renders  A  and  B  in  (2.18)  diagonal 
simultaneously,  the  matrix  L(T;d)  will  become  diagonal  and  the  HCM  will  converge  in 
one  step. 

Corollaiv  4.3.  Let  A  =  | '^b(z)^\^^dz,  B  =  | with  the  Junctions  yj/ 
defined  in  (2.13)-(2.15).  Let  Q  be  the  (q+1)  x  (q+1)  matrix  the  coiumns  of  which  are  the 
eigenvectors  of 

(^•27)  Bx  =  afAx 

normalized  so  that  x^Ax  =  1  and  0  ^  a?  ^  of,,.  Then  the  basis  ip  =  renders 
L(ijd)  in  (4.14)  diagonal  and  the  HCM  converges  in  one  step,  ic.  x  =  0  in  (3.5). 

Proof.  We  have  x^x  =  x'^Ax  =  <rf  by  (4.19).  Now,  due  to 

|’  b(z)^V»^dz  =  Q^AQ  «  1, 

J^^z)y»V'^dz  =  Q^BQ  =  {0,aj,...,a;}, 

hence  with  =  0 

L(f)  •  f=A  ♦  B  »  diag(f'*i.?}’_^.  □ 

Frequently,  especially  in  elasticity  problems,  it  is  not  possible  to  find  a 
independent  transformation  Q  which  renders  L  diagonal.  Nevertheless,  as  we  shall  see  in 
the  next  section,  it  is  often  useful  (and  possible)  to  have  L(^)  approach  diagonal  form 
as  f  -*  0*  or  as  f 

Coiollarv  4.4.  Assume  that  0  <  ^  <  2,  Then 
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V:  lim  \(6,ri)  »  0  if  B  =  f  ’  a(z)ip' (ip' y  dz  is  diagonal 

ij-*0  J 

2®:  lim  =0  if  A  =  f*  h(7) tpt^dz  is  diagonal 

Tf^CO  J 

Proof.  1®.  Consider  L(f)  as  f -*  0.  Since  then  +  B  -►  B  and  B  is 

diagonal,  the  assertion  would  follow  if  B  were  nonsingular.  Since,  however,  rank  (B)  = 
q  <  dim(B)  =  q  +  1,  we  assume  that  ol  =  0,  <p^  =  const,  and  that 
’^b(z)^j^Qdz  B  0  for  i  2:  1.  This  uncouples  the  first  equation  in  (4.18)  from  the 
remaining  ones,  i.e.  the  iteration  matrix  attains  block  diagonal  form: 

= 

The  corresponding  lower  right  block  of  B  is  now  diagonal  and  positive  definite,  hence 
1®  follows. 

2®.  We  write  in  this  case 

UD  ‘  «A-r‘B). 

From  the  formula  for  M  in  Theorem  4.1,  2®  we  find  that  M  corresponding  to  L  in 
(4.20)  and  M  corresponding  to  (A  j^^B)  are  identical  (i.e.  a  scaling  of  L  does  not 
alter  M).  Since  A  is  diagonal  and  positive  definite,  and  f"^B  -»  0  as  f  -♦  oo,  we 
find  2®.  □ 

In  summary,  we  have  shown  that  the  influence  of  the  basis  on  the 


25 


convergence  rate  of  the  HCM  Algorithms  3.1  can  be  analyzed  quaiititatively  with  the 
iteration  matrix  in  (4.15).  In  the  same  vein  any  other  iterative  method  that 

is  based  on  the  subspace  decomposition  (3.1)  can  be  analyzed,  in  particular  the  parallel 
subspace  correction  method. 

Remark  4.4.  The  analysis  presented  so  far  applied  only  to  the  case  where  a>  is  an 
interval  or  a  rectangle  and  the  boundary  conditions  are  either  homogeneous  Dirichlet-  or 
periodic  conditions.  However,  analogous  results  can  be  obtained  for  arbitrary  Lipschitz- 
domains  «  C  R”,  if  the  unknown  coefficient  functions  X(x)  and  their  iterates  X®(x) 
are  expanded  into  series  of  eigenfunctions  ^^(x)  of  the  operator  -A^  in  <■>  with 
appropriate  homogeneous  boundary  conditions  on  du>. 

S.  FOURIER  CONVERGENCE  ANALYSIS  OF  THE  HCM-U 

In  the  previous  section  we  introduced  the  tool  of  Fourier  transformation  of  the 
reduced  model  (3.3)  to  obtain  explicit  estimates  on  the  convergence  rate  k  in  (3.5)  and 
p  in  (3.6).  This  was  achieved  by  taking  the  supremum  overall  wavenumbers  m  of  the 
largest  eigenvalue  X(0,ij„d)  of  (4.16),  resp.  of  r(M(^,i;_d))  of  (4.17).  Our  purpose  in 
the  present  section  is  to  show  that  the  eigenvalue  curve  resp.  r(M(0,Yr)) 

contains  more  information:  it  is  a  factor  by  which  the  energy  contained  in  error 
components  of  "wavelength"  i;'‘d  is  reduced  in  one  step  of  the  HCM.  We  consider 
here  the  special  case  of  (2.5)  where 
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(5.1) 


a  =  b  -  1,  C  =  1  (unit  matrix), 


so  that  L  =  A 
in  (2.13)-(2.15) 


a 

and  =  —  Then  it  is  readily  verified  that  the  director 
dn 

are  even  polynomials  of  z  and  we  list  the  first  of  them  in 


functions  rp. 
the  following 


table. 


T/BLE5.1.  The  function  ^^(z),  0  <  j  !£2. 


j 

^j(z) 

0 

1 

1 

(3z"  - 1)/6 

2 

(15z"  -  30z^  +  z)/360 

The  basis  functions  which  render  A  in  (2.17)  diagonal  are  orthogonal  in 
L*(-l,l)  and  we  select  them  to  be  the  even  Legendre  polynomials,  i.e. 

(5.2)  ^z)  .  L,(z) 

where  Lo(z)  =  1,  L2(z)  =  (3z*-l)/6,  etc. 

Analogously,  the  functions  which  render  B  in  (2.17)  diagonal  are  orthogonal  in 
H’(-l,l)  with  respect  to  the  inner  product 

(5.3)  =  |‘^'./.dz. 

A  basis  which  satisfies  (^j,^j)  =  is  therefore  given  by 
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(5.4) 


^o(2)  =  h  V>j(2)  = 


±1 

2 


or,  more  explicitly,  <p.(z)  «  (2(4j  -  l))->«(L5j(z)  -  Lj..j(z)). 

Figures  5.1  and  5.2  depict  the  bounds  p(M($,ti))  for  0  -  1  (i.e.  no  relaxation)  over 
several  orders  of  magnitude  of  the  normalized  frequency  i;  for  the  basis  function  in  (5.2) 
and  (5.4),  respectively.  We  remark  that  Lo(2)  and  Lliz)  are  orthogonal  in  the  7  C-l|l) 
and  in  the  H*(-l,l)  inner  product  (5.3)  -  hence  we  have  X(6,ij)  =0  for  q  =  1  with  the 
basis  (5.2),  in  accordance  with  Corollaiy  4.3. 

In  accordance  with  Corollaiy  4.4,  we  have 

(53)  lim  p(M(l,i,))  =  0  for  (5.2) 

and 

(5.6)  lira  p(M(l,i,))  =  0  for  (5.4), 

i.e.  (5.2)  is  better  for  a  rapid  convergence  on  the  high  frequency  components  while  (5.4) 
performs  superior  on  the  low  frequencies.  Of  principal  interest  is  the  overall  supremum 
for  ®*1,  q  =  1, 2,3,4  with  the  basis  (5.4). 


(5.7)  sup  (p(M(1,7))}  p(q). 

In  Table  5.2  we  report  the  numerical  values  for  the  bases  (5.2)  and  (5.4). 
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TABLE  5.2:  p(q)  in  (5.7)  for  the  bases  (5.2)  and  (5.4). 


q 

0 

1 

2 

3 

4 

(5.2) 

0 

0 

0.3000 

0.4764 

0.5833 

(5.4) 

0 

0.833 

0.9400 

0.9694 

0.9811 

We  see  clearly  that  the  basis  (5.2)  is  more  robust  than  (5.4),  since  it  yields  reasonably  small 
contraction  rates  uniformly  on  all  frequencies  The  performance  of  the  basis  (5.4)  can 
be  poor  if  the  error  contains  high  frequency  components.  In  the  following  table  we  list 
for  q  =  1  and  the  basis  (5.4),  for  small  i). 


TABLE  5  J.  p(q)  for  q  =  1  and  the  basis  (5.4). 


V 

0.001 

0.01 

0.1 

0.25 

0.5 

p(M(1,i7)) 

0.33*  10-* 

0.33- 10" 

0.22.10-' 

0.7510' 

We  see  clearly  that  we  have  here 

(^•®)  p(M(l,»j))  ~  ij^/3  as  ij  -♦  0*. 

In  Figure  5.3  we  depict  this  quantity  also  for  higher  q,  with  the  slope  according  to 

(5.8). 
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Figure  5.3.  p(M(tf,ij))  for  e  =  1,  q  =  1,23,4  with  basis  (5.4),  log-log  scaie. 

V  ce  further  that  (5.8)  is  practically  independent  of  the  model  order  for  ii  ^  1. 

The  above  graphs  show  that  we  can  expect  generally  a  different  performance  of  the 
HCM  on  high  and  low  frequency  components  of  the  error  Y^^x)  :=  X(x)  -  X^''^*).  This  is 
of  significance  since  the  exact  solution  X(x)  consists,  especially  for  small  d,  of  two 
mutually  distinct  components:  a  smooth,  asymptotic  part  and  houmfiuy  ^^(x)  of  the 
form 

(5-9)  X,“(x)  . 

where  depends  only  on  the  model  order  q  [16]  and  xf^(x)  is  smooth. 
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In  Table  5.4  we  list  the  values  of  oj(q). 


TABLE  5.4;  Boundary  layer  «q)onents  (rj(q)  for  small  q. 


q 

1 

2 

3 

oo 

‘r,(q) 

3.1529 

3.1416 

T 

‘r2(q) 

9.7498 

6.4791 

2t 

‘TaCq) 

18.0596 

3t 

It  is  further  shown  in  [16]  that 

ojCq)  -*  j*"  as  q 

Comparing  formally  exp(-<rx/d)  and  exp([m/3]x)  with  jS  =  -  we  expect  large  Fourier 

2 

coefficients  x,  for  [m/J]  -  <r/d  or  for 


(S.10)  V  - 

a 

We  remark  at  this  point  that  for  certain  compatible  data  f(x)  the  components  (5.9) 
vanish.  Nevertheless,  for  arbitraiy  data  the  boundaiy  layers  (4.9)  are  generally  present  and 
we  see  from  (5.10)  that  we  can  expect  the  contraction  rate 


(5.11) 


P  ^  X 


X 

G,—] 

«J 

for  the  HCM.  If  several  distinct  layers  are  present  (cf.  Table  5.4),  we  see  that 
or/a  can  be  quite  large,  hence,  a  decay  of  X(9,i;)  as  ri  -*  oo  is  highly  desirable. 
Remark  5.1.  It  follows  from  Theorem  4.1  only  that  X  ^  k.  Nevertheless,  we  found  in  all 


cases  for  the  matrix  M  in  (4.15)  that  p  «  X. 


6. 


NUMERICAL  EXPERIMENTS 


In  order  to  verify  the  applicability  of  our  theoiy  in  a  computational  setting,  we 
considered  the  following  model  problem: 

0  -Au  *  0  in  (-1,1)  X  (-d,d), 

(^•^)  u  =  0  for  X  =  ±1 

(6.3)  ^  =  ±1  for  y  =  ±d. 

dy 

The  solutions  exhibit  boundaiy  layers  of  the  form  (5.9)  (see  also  [16]). 

The  problem  (6.1)-(6.3)  was  discretized  by  a  two  dimensional  FBM  with  rectangular 
elements  Wj  x  (-d,d)  where  Uj  C  (-1,1)  are  subinteivals.  Tensor  product  polynomial 
subspaces  with  degree  p  in  horizontal  and  uniform  degree  q  (only  even  shape  functions 
were  used  here)  in  the  transverse  direction  were  used.  We  selected  p  s  10  to  ensure  that 
the  models  were  solved  accurately. 

Remark  6.1.  In  practice  a  locally  variable  degree  q  is  of  high  interest,  see  e.g.  [13].  In  this 
case  the  HCM  becomes  a  block  iterative  method  (here  a  block  SOR)  for  the  fully  discrete 
system. 

An  exact  solution  was  obtained  by  direct  solution  of  the  linear  ^tem  and  the 
iteration  was  stopped  when 


(6.4) 


<  10-* 

I 


where  ||  •  ||  e  denotes  the  (discrete)  energy  norm.  We  selected  the  initial  guess  =  0 
throughout  (in  practice  when  solving  the  higher  order  models,  one  should  start  from  a 
converged  solution  of  the  lower  order  problem).  As  an  averaged  convergence  rate  we  took 
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the  geometric  mean  of  J  iterates  before  convergence  occurred,  i.e. 


(6J) 


f; 

k-J-,  ^-‘'-x^IIe 


i/j 


(J  =  min{10,k'}  was  selected).  In  Figure  6.1  we  show  the  average  convergence  rates  for 
models  based  on  the  director  functions  (S.2).  The  model  orders  q  =  0  and  q  =  1  were 
not  shown,  since  the  HCM  converges  in  one  step  here. 

We  observe  in  Figure  6.1  that  the  maximum  convergence  rates  are  veiy  close  to  the 
values  predicted  by  our  Fourier  analysis.  Moreover,  we  observe  that  as  d  0  all  rates 
tend  to  0.3,  the  theoretical  value  (cf.  Fig.  5.1).  This  is  due  to  the  contribution  of  the  q  = 
3  and  4  models  being  of  higher  order  as  d  0  and  due  to  the  1-step  convergence  of 
HCM  for  q  =  1. 

Let  us  now  turn  to  the  HCM  with  the  basis  (5.4).  Figure  6.2  shows  the  performance 
for  (6.1)-(6J). 

Note  that  for  the  boundary  layer  case  (Fig.  6.2),  the  convergence  rates  are  once  again 
close  to  the  ones  predicted  in  Table  5.2  which  are  in  this  case,  for  large  d,  very  close  to 
1. 

Finally  we  present  in  Figure  6.3  the  convergence  history  for  the  HCM  with  (5.2)  for 
various  values  of  d  and  for  q  =  2.  Figure  6.4  presents  the  analogous  data  for  the  basis 
(5.4). 
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RELATIVE  ENERGY  NORM  ERROR 


NUMBER  OF  ITERATIONS 


7. 


CONCLUSIONS 


We  have  shown  that  the  hierarchic  models  can  be  solved  very  effectively  by  an 
iterative  procedure,  provided  that  the  shape  functions  used  (in  y  direction)  are  properly 
selected.  The  iterative  procedure  employed  was  based  on  a  subspace  decomposition  with 
respect  to  the  spectral  order  of  the  hierarchic  plate  model.  Shaip  estimates  for  the 
convergence  rate  of  the  subspace  correction  in  dependence  on  the  selection  of  the  director 
functions  were  obtained  using  the  Fourier  transform.  The  technique  will  be  extended  to 
elasticity  problems  in  the  forthcoming  second  part  of  this  work. 
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